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Abstract 

We consider FPU-type atomic chains with general convex potentials. The naive continuum hmit 
in the hyperbolic space-time scaling is the p-system of mass and momentum conservation. We sys- 
tematically compare Riemann solutions to the p-system with numerical solutions to discrete Riemann 
problems in FPU chains, and argue that the latter can be described by modified p-system Riemann 
solvers. We allow the flux to have a turning point, and observe a third type of elementary wave 
(conservative shocks) in the atomistic simulations. These waves are heteroclinic travelling waves and 
correspond to non-classical, undercompressive shocks of the p-system. We analyse such shocks for 
fluxes with one or more turning points. 

Depending on the convexity properties of the flux we propose FPU-Riemann solvers. Our nu- 
merical simulations confirm that Lax-shocks are replaced by so called dispersive shocks. For convex- 
concave flux we provide numerical evidence that convex FPU chains follow the p-system in generating 
conservative shocks that are supersonic. For concave-convex flux, however, the conservative shocks 
of the p-system are subsonic and do not appear in FPU-Riemann solutions. 

1 Introduction 

The derivation of effective continuum descriptions for high-dimensional discrete systems is a fundamental 
tool for model reduction in the sciences. Hamiltonian lattices, such as atomic chains, naturally lead 
to nonlinear systems of conservation laws which describe the leading order dynamics on the hyperbolic 
space-time scale. It is customary to neglect the higher order terms and to study the leading order system 
by itself. The rigorous mathematical validity of this reduction is a notoriously difficult task and there 
are surprisingly few successes reported in the literature. 

This paper concerns the macroscopic description of monoatomic chains with nearest neighbour inter- 
actions, see Figure [TTTI In their seminal paper [FPU55| Fermi, Pasta and Ulam chains studied such chains 
for interaction potential whose non-harmonic part involves only cubic or quartic terms. We allow for 
general convex interaction potentials, but still refer to the systems as FPU chains. 
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Figure 1.1: The atomic chain with nearest neighbour interaction. 



An important building block for the macroscopic descriptions of FPU chains are the solutions to 
Riemann initial data in the hyperbolic continuum limit. In this limit space and time are scaled in the 
same way, and the amplitude of solutions is unconstrained. Since the numerical study of Holian and 
Straub |HS78| and from rigorous results for the integrable Toda chain it is known that the solutions to 
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atomistic Riemann problems with either convex or concave flux obey a self-similar structure on the 
hyperbolic scale: each solution consists of at most two elementary waves that are separated by constant 
states. Each of these elementary wave is cither a rarefaction wave, in which the atomic data vary smoothly 
on the macroscopic scale, or a dispersive shocks, in which strong microscopic oscillations spread out in 
space and time. 

The starting point for our investigation was the observation that for certain $ a third kind of elemen- 
tary waves can be observed in FPU-Riemann problems. These waves, which we refer to as conservative 
shocks, involve no oscillations and look like 'shocks' (jump discontinuities) on the macroscopic scale. Of 
course, these waves are not exact shocks as they exhibit a transition layer on the atomistic scale, but 
this layer is very small and disappears in the hyperbolic scaling. To our knowledge the appearance of 
conservative shocks in FPU-Riemann problems was never reported before. 

An illustrative example of an FPU-Riemann problem is plotted in Figure 11.21 and involves all types 
of elementary waves (from left to right): a rarefaction wave, a dispersive shock, and conservative shock. 
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Figure 1.2: Riemann problems in convex FPU chains can involve three kinds of elementary waves: rar- 
efaction waves, dispersive shocks, and conservative shocks. The first two pictures show snapshots of the 
atomic distances and velocities against the scaled particle index ajN for a chain with N ~ 8000 particles 
and flux $'(r) = ^(r ^ arctanr); the third pictures magnifies the region around the conservative shock. 

The naive approach to describe the FPU dynamics on the hyperbolic scale assumes long-wave-length 
motion without microscopic oscillations. Under this assumption one readily derives the 'p-systcm', which 
consists of the conservation laws for mass and momentum in Lagrangian coordinates. 

It is well known, that the naive continuum limits of nonlinear dispersive lattices provide a reasonable 
macroscopic model as long as the macroscopic fields are smooth, see, e.g., |Lax86[ [GL88[ IHL9T1 IHLM94[ 
ILL96| . The nonlinearity, however, usually causes shock phenomena, and the naive continuum limit fails 
in this case. Instead, lattice systems like FPU typically produce dispersive shocks, in which the atoms 
self-organize into strong microscopic oscillations. On the macroscopic scale such dispersive shocks can be 
regarded as measure- valued solutions to the naive continuum limit, see (|2.9|) in fj2j 

The formation of dispersive shocks is a characteristic property of Hamiltonian 'zero dispersion' limits 
and a direct consequence of the conservation of energy, see ^ Moreover, it is known from numerical stud- 
ies and rigorous results for integrable systems, see |GP73[ |LL83[ IVen85l ILax86[ ILaxQll ILLV93[ IKamOO] . 
that the oscillations in a dispersive shocks arc modulated wave trains (period travelling waves). Figurc fLBl 
presents a typical example of a dispersive shock in FPU chains. At the shock front, where the amplitudes 
of the oscillations become maximal, the wave trains converge to a supersonic soliton, that is a homoclinic 
travelling wave. 

By combining non-classical hyperbolic theory of the p-system with macroscopic theory, travelling 
waves and numerical observations of FPU chains we characterise FPU-Riemann solvers for oscillation- 
free initial data and fiuxes with one turning point. 



Conservative shocks in FPU chains and the p-system. As shown in Figure [TT^ the solution 
to FPU-Riemann problems with convex $ can involve conservative shocks. Below in ^we explain that 
these waves correspond to certain shocks in the p-system, namely those that conserve the energy exactly. 
Among all p-systcm shocks, the set of conservative shocks is quite small and if has no turning point 
conservative shocks cannot occur at all. The conservative shocks in the p-system are non-classical shocks 
as they violate the Lax-condition: for convex-concave there are fast undercompressive, and hence 
supersonic with respect to both the left and right state; for concave-convex $' the conservative shocks 
are slow undercompressive and hence subsonic. 
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Figure 1.3: Dispersive shocks arise naturally in FPU-Ricmann problems as a consequence of energy 
conservation. Each dispersive shock is build of a one-parameter family of wave trains with a single soliton 
at the leading front. (a,b) Snapshots of atomic distances and velocities, and their local mean values, (c) 
Superposition of several local distribution functions within the shock; positions of the mesoscopic space- 
time windows are marked by vertical lines in (a,b). 



We numerically discovered that conservative shocks occur naturally in FPU-Riemann problems chains 
if they are supersonic. However, numerical simulations with concave-convex $' never generated anything 
close to a jump discontinuity. Instead, we typically observe solutions as plotted in Figure [L^ the solution 
appears to be a composite wave of a dispersive shock with attached rarefaction wave. 

Conservative shocks in FPU chains are naturally related to atomistic fronts which are 'heteroclinic' 
travelling waves. A bifurcation result of looss [looOOj shows the existence of small amplitude atomistic 
fronts for convex-concave flux (in which case the conservative shocks are supersonic) . The authors have 
improved this result in [HR09| by showing that each front must correspond to a conservative shock in 
the p-systcm and that supersonic shocks with arbitrary large jump height can be realised by an atomistic 
front. 
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Figure 1.4: FPU-Ricmann problem with initial data corresponding to a subsonic conservative shock in the 
p-system: Instead of a conservative shock the atoms self-organize into a dispersive shock with attached 
rarefaction wave. The leading soliton in the dispersive shock is hence sonic, (a) atomic distances; (b) 
atomic velocities; (c) five local distributions hmctions corresponding to the vertical lines in the snapshots. 



FPU-Riemann solvers. A main goal of this paper is the derivation of FPU-Riemann solvers which 
predict the number and the type of the elementary waves that result from arbitrary Riemann initial 
data. To this end we systematically compare numerical simulation for various FPU chains with certain 
Riemann solvers for the p-system. The FPU-Riemann solver for the Toda chain is well understood, see 
for instance |Kam93[ IDM98| , but there is no complete picture for non-integrable chains as the available 
numerical studies only concern special types of Riemann initial data, such as the 'piston problem' in 
[HS78j . Moreover, we are not aware of any previous analytical or numerical investigation of any FPU- 
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Riemann problem which allows for conservative shocks. 

In the classical case with cither convex or concave flux the numerical simulations indicate that 
the solution to each FPU-Riemann problem can be described by an adapted classical solver, in which 
Lax-shocks are replaced by dispersive shocks. More precisely, from each given left state there emanate 
four curves (wave sets) in the state space. Two of them correspond to rarefaction waves and appear 
also in each solver for the p-system. Instead of Lax-shock curves, however, we find two dispersive shock 
curves. The solution to each FPU-Riemann problem is then completely determined by these wave sets. 
In particular, in the classical case each Riemann solution consists (generically) of a left moving 1-wave 
and a right moving 2-wave, where each wave is either a rarefaction wave or a dispersive shock. 

In presence of turning points of <&' the Riemann solvers for both FPU chains and the p-system are 
more complicated because now the solutions to general Riemann problems involve composite waves, which 
consist of two elementary waves from the same family, and may also involve undercompressive shocks. The 
Riemann solvers in these cases can be described in terms of modified wave sets, but for the p-system they 
are not unique. Consequently, different p-system Riemann solvers are possible, such as the 'conservative' 
and the 'dissipative' solver described in SjH For FPU chains, however, the underlying atomistic dynamics 
determines the solver uniquely for each 




Figure 1.5: Sketch of the modified rarefaction curve (rp{ < tl) for convex-concave $' and given tl, and 
corresponding FPU-Riemann solutions: For < tl one still finds a rarefaction wave (R), but if tl 
crosses the turning point of $' a non-classical, a conservative shock (N) nucleates. If tr decreases 
further the rarefaction wave is replaced by a dispersive shock (D). 

In we investigate systematically numerical solutions to FPU-Riemann problems with convex- 
concave and concave-convex $' and argue that the solutions can be described by adapted p-system 
solvers. More precisely, in the convex-concave we propose to adapt the conservative solver, which pre- 
dicts composite waves with supersonic conservative shocks. The modified rarefaction curve of such a 
solver is illustrated in Figure 11.51 In the concave-convex case, however, the simulations indicate that 
FPU-Riemann solutions can be described by an adapted dissipative solver, whose modified wave sets 
do not involve conservative shocks. In both non-classical FPU solvers the Lax shocks are replaced by 
dispersive shocks, and in the nucleation criterion for composite waves the dispersive shock front velocity 
plays the role of the Rankine-Hugeniot condition. These adaption rules lead to an adequate description of 
the numerical solution to FPU-Riemann problems. Moreover, the difference between Rankine-Hugeniot 
and dispersive shock front velocity provides an explanation for the absence of conservation shock in the 
concave-convex case: the nucleation criterion for conservative shocks cannot be satisfied. However, it is 
not known how to predict the shock front velocity from the initial data. 

We believe that these results give new structural insight into the hyperbolic nature of the continuum 
limit and the role of the p-system in it. Moreover, they open up new avenues for analytical investigation, 
some of which we phrase in the form of conjectures. 

We emphasize that the situation dramatically changes for fluxes with more than one turning point. On 
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the p-system level, non-classical Riemann solvers can still be constructed, but these no longer provide a 
basis for the continuum limit of the FPU chain. The reason is that we numerically find energy conserving 
shocks between wave trains. This is a new phenomenon but its study is beyond the scope of this article. 
Another problem left for future investigation concerns Riemann problems where oscillations in form of 
wave trains are already imposed in the initial data. One of the expected new phenomena in this case is 
the onset of two-phase oscillations. Finally, it would be interesting to study cold initial data with more 
than one jump discontinuity, because then two-phase oscillations can be created by the interaction of two 
dispersive shocks. 

Organisation of the paper. In fj^l we collect some facts about FPU chains and the p-system; 
we especially discuss how oscillatory FPU data give rise to measure-valued solutions to the p-systems 
and briefly outline the concept of modulated wave trains. 331 concerns the numerical simulation of 
FPU-Riemann problems and contains our observations about self-similarity on the macroscopic scale; 
we also investigate the fine-structure of dispersive and conservative shocks. §4] is devoted to FPU- 
Riemann problems. We start with numerical results for the classical case with either convex or concave 
Afterwards we discuss the conservative and dissipative solvers for the p-system, and proceed with 
FPU-Riemann problems in the non-classical cases. Finally, in fJS]we prove some analytical results about 
conservative shocks in the p-system. 

2 Preliminaries of FPU chains and the p-system 

Convex FPU chains consist of N identical particles which are nearest neighbour coupled in a convex 
potential $ : K ^ M by Newton's equations 

Xa = ^'{Xa+l - Xa) - ^'{x^ - Xa-l), (2.1) 

where ' = ^ is the time derivative, Xa{t) the atomic position, a = 1,...,N the particle index. We 
consider nonlinear force referred to as the flux. A prominent example for nonlinear (without 
turning points) is the completely integrablc Toda chain, see |Tod70[ IFla741 IHen74| , with 

$Toda(r) =exp(l-r)-(l-r). (2.2) 

For our purposes it is convenient to use the atomic distances = and velocities 

the basic variables, changing (|2.ip to the system 

Ta = Va+l - Va , = $' (r^ ) - $' (ra_ i ) . (2.3) 

Note that while the Toda potential, and also the potentials used below, allows for negative distances, this 
is essentially a matter of suitably shifting the minimum by $ ^ $(• 4- tq). 

We are interested in the thermodynamic limit £ = l/iV^Oin the hyperbolic scaling of the microscopic 
coordinates t and a. This scaling is defined by the macroscopic time t = et and particle index a = ea. 
It is natural to scale the atomic positions in the same way, i.e. x = ex, which leaves atomic distances 
and velocities scale invariant. In the limit e = the spatial variable a becomes continuous and the high 
dimensional ODE p.ip should be replaced by a continuum limit, i.e., by a system of a few macroscopic 
PDEs. Microscopic oscillations can be naturally interpreted as a form of temperature in the chain, see 
jDHR06| , and accordingly we refer to oscillation-free limits as cold. 

Evolution of cold data and the p-system To derive the p-system as the simplest model for 
the macroscopic evolution of FPU chains we assume macroscopic fields r(i,a) and v{t,a) such that 
Tait) = r{et, ea), Va{t) = v{et, ea). This ansatz corresponds to cold motion as it assumes that there are 
no microscopic oscillations in the chain. Substitution into (|2.3p and taking the limit e ^ yields the 
macroscopic conservation laws for mass and momentum 

djr-doiv = 0, djv - da<S^'{r) ^0. (2.4) 
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It is well known that the p-system is hyperbolic for convex $ and that for smooth solutions the energy 
is conserved via 



dj + $(r)) - da {v $'(r)) = 0. 



In the p-system a shock propagates with a constant shock speed Crh so that r and v satisfy the Rankine- 
Hugeniot jump conditions for mass and momentum 



where \[x]\ ^ xi^ — denotes the jump. The main observation is that for either convex or concave flux 
<&' the jump conditions (|2.5p imply that the jump condition for the energy must be violated, i.e., 



see Theorem 15.11 below. Here the Lax criterion selects the shocks with negative production. 

Onset of dispersive shocks in FPU chains It is known that the p-system provides a reasonable 
thermodynamic limit for FPU chains in the following sense: Preparing cold FPU initial data with smooth 
profile functions r and v, the atomistic dynamics reproduces a solution to the p-system provided that the 
latter has a smooth solution. This can be understood as a manifestation of Strang's Theorem [Str64| : we 
refer to |DH08| for numerical simulations and to [GL88. ,HL91, for a similar discussion in the context of 
other lattice systems. 

At some critical time, however, the p-system forms a shock; this shock still conserves mass and 
momentum according to (|2.5p . but in most cases it has a negative energy production. In contrast, 
Newton's equations always conserve mass, momentum and energy, so the continuum limit of FPU chains 
beyond a shock cannot be described in terms of the p-system. Instead, the FPU chain produces a 
dispersive shock with strong microscopic oscillations that take the form of modulated wave trains, compare 
Observations 13.21 and 13.31 below. 

Heuristically, the onset of dispersive shocks is a consequence of energy conservation and can be 
interpreted as Hamiltonian self-thermalisation: When the shock is formed the p-system predicts some 
macroscopic excess energy which no longer can be stored in cold motion. On the atomistic scale this 
excess energy is transferred into modulated wave trains and appears as internal or thermal energy on 
the macroscopic scale. More precisely, although wave trains do not provide a thcrmalization in the usual 
'chaotic' sense, their macroscopic dynamics is governed by thermodynamically consistent field equations, 
see |DHM06llDHR06] . 

Riemann solvers for the p-system The p-system is hyperbolic if $ is convex and genuinely 
nonlinear for 7^ 0, thus turning points of $' correspond to states in which the system is linearly 
degenerate. The Riemann problem for strictly convex or concave flux can therefore be described by 
the classical solver, which is based on the classical Lax theory from |Lax57| and involves only rarefaction 
wave and Lax shocks. 

This classical Riemann solver is built from the following curves, where — corresponds to left moving 
1-waves and + to right moving 2-waves. The rarefaction wave sets TZ± [ml] contain all right states 
UR = (^R7 vb.) that can be reached from a given left state wl ~ i'l) with a single 1- or 2-rarcfaction 
wave. The shock wave sets S±[ui] consist of all possible right states wr that can be reached by a single 
Lax 1- or 2-shock. The sets yV-|-[uL] = 7?.±[ul] U 5±[ul] form C^-smooth curves through ul, and we 
denote W [wl] = W- [ul] U yV+ [ml] . The solution to the Riemann problem with given left and right states 
Ml and mr consists of the two elementary waves that connect itL to mm, and um to mr (one of these may 
be trivial), where the intermediate state is uniquely determined by mm G W_[ml] and mr G yV-)-[MM]- 

In Appendix |^ we give more details about the classical-solver for the p-systcm. If $' has turning 
points the wave sets of the classical solver must be modified, and this gives rise to non-classical solvers 
which involve various types of composite waves, see m.2\ 

Conservative shocks in the p-system In contrast to Lax shocks, conservative shocks in the 
p-system balance mass, momentum p.Sp and energy (|2.6p . This gives rise to the system of nonlinear 



cM + \[v]\^0, cM + Wir)]\=0, 



(2.5) 



c,4y + ^r)]\ + \[v^'{r)]\^0, 



(2.6) 
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equations 



crh W + HI = 0, cM + W{r)]\ = 0, c,4h' + '^(^)]\ + = 0' (2-7) 

for the five parameters tl, wl, ''r, t'R, Crh- According to Appendix [Xl a conservative shock is called 
supersonic if |crh| > inax{|A±(rL)|, |A±(rp{)|} and subsonic if |crh| < inin{|A±(rL)|, |A±(rR)|}. It can be 
easily shown that each conservative shock satisfies 

J(rL, rR.) := |[a>(r)]| - |[r]| (<i>'(r)) = 0, ($'(r)) i($'(rL) + <i>'(rH)). (2.8) 

Conversely, for each solution to (|2.8p there exist both a corresponding conservative 1-shock and 2-shock. 
These shocks are unique up to Galilean transformations, and differ only in sgn|w]| = sgncrh- We analyse 
the set of conservative shocks in the p-system in more detail in ^ 

Macroscopic limit and Young measures About the hyperbolic continuum limit of FPU chains 
in the presence of strong microscopic oscillations little is known rigorously. This is the reason why, 
for a large part, we have to rely on numerical observations. The main difficulty lies in the control of 
oscillations that lead to measure-valued solutions on the macroscopic scale. Heuristically, such measure- 
valued solutions are governed by extended p-systems, but a rigorous derivation of such extensions could 
be treated in a satisfactory manner only for integrable systems so far; notably the harmonic chain 
|DHM06[ IMieOel IMac02[ IMac04| . the hard sphere model [HcrOSj, and the Toda chain |DKKZ96[ [DM98] . 

Nevertheless, some insight into the macroscopic evolution of microscopic oscillations can be gained 
from the theory of Young measures. Here it is supposed that the atomic data generate a family of 
probability distributions fi(t, a, dQ) , where (t, a) is a point in the macroscopic space-time and Q = (r, v) 
denotes a point in the microscopic phase space of distances and velocities. Note that the atomic data 
are oscillation free in the vicinity of a point (i, a) if and only if the measure ii(t, a, dQ) is a delta 
distribution with respect to the Q variable. 

On the one hand. Young measures provide an elegant framework to investigate oscillatory numerical 
data that we used to interpret our simulations. For given (t, a) the measure a, dQ) can be approxi- 
mated by means of mesoscopic space-time windows and provides local mean values of atomic observables 
as well as statistical information about the microscopic oscillations, see |DHM06[ IDH08j . 

On the other hand, Young measures are useful for analytical considerations because the solutions to 
(|2.ip with iV — > oo are compact in the sense of Young-measures provided that the initial data are of order 
1. Extracting convergent subsequences, one can then prove as in [HerOSj that every limit measure must 
be a weak solution to the following macroscopic conservation laws of mass, momentum and energy 

dj{r)^dr,{v)^0, 
dj{v)-d^{^'{r))^0, (2.9) 
dj{^v^ + $(r)) - dr, {v $'(r)) ^ 0. 

Here (4')(i, a) is the local mean value of the observable '5 = 'I'(r, v), that is 

(*)(t,a) = / lPiQ)^l{la,dQ). 

System (|2.9p provides non-trivial information about the macroscopic dynamics of FPU chains. For 
arbitrary oscillations, however, we can not express the fluxes in terms of the densities, and hence (j2.9p 
does not determine the macroscopic evolution completely. 

As an important consequence of (|2.9p we can characterise conservative shocks in FPU chains. To this 
end suppose that for all points (i, a) in a sufficiently small region the measure /i(i, a, dQ) depends only 
on c = a/t and is a delta-distribution Sui^{dQ) for c < and 6u^{dQ) for c > Crh, for some Crh- This 
is exactly what wc observe in the numerical FPU simulations for a conservative shock connecting ul to 
Mr, see Figure [L2l In this case ()2.9p reduces to the three independent jump conditions that determine a 
conservative shock in the p-system, compare p.7p . In other words, FPU chains allow for waves that are 
close to a jump discontinuity only if there exists a corresponding conservative shock in the p-system. 
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Travelling waves and modulation theory Careful investigations of numerical experiments as 
described in [DH08j reveal that also for non-integrable cases the oscillations in a dispersive shock of FPU 
chains can be described by modulated travelling waves. That means for each (t, a) in the oscillatory 
region the measure a, dQ) is generated by a travelling wave, whose parameters are slowly varying as 
they depend only on the macroscopic coordinates t and a. This observation is in accordance with the fact 
that the support of /i(<, a, dQ) is contained in a closed curve, compare the density plots in Figures 11.41 
and ll.51 which show the superposition of several of these curves. 

Travelling waves with constant speed c arc exact solutions to the infinite chain (|2.ip that depend on 
a single phase variable = ka + cot via Xa(t) = x{(j)). Here k and lo arc generalized wave number and 
frequency, respectively, and c = —uj/k is the phase velocity. In terms of atomic distances and velocities 
travelling waves can be written as ra{t) = R{4>) and Va(t) = V{(f>), where the profile functions R and V 
solve the advance-delay differential equations 

ca^i?(0) = F(0 + 1) - Vi^) c9^F(0) ^ $'(i?(0)) - ^'{Ri^ - 1)). 

In our context relevant travelling waves are wave trains, for which both R and V arc periodic functions, 
solitons (solitary waves), which limit to the same background state as — > ±cx), and fronts, which connect 
to different constant background states for (/> — > — oo and </) — > oo. 

Wave trains exist for all convex potentials 4>, see jFV99[ IDHM06I IHerOSj . They depend on four 
parameters and provide the building blocks for modulation theory, which describes the macroscopic 
evolution of a modulated travelling wave. This evolution is governed by a system of four nonlinear 
conservation laws, which one usually refers to as Whitham's modulation equations, and a dispersive 
shock is just a rarefaction wave of this system. The Whitham equations can be regarded as an extension 
of (j2.9p . where modulation theory also provides a complete set of constitutive relations which depend 
on $ via the four parameter-family of wave trains. We refer to [Whi74| for the general background, and 
jFV991 IDHM06j for the apphcation to FPU chains. 

The existence of solitons for super-quadratic potentials is proven in broad generality in jFW94|[SW97| . 
and [PPOOI IHerOSj show that wave trains limit to solitons as the wave number tends to zero, see also 
[PanOSj . Solitons are important in our context as they appear at the shock front of a dispersive shock 
where the amplitude of the oscillations is maximal. Generically solitons travel faster than the sound 
speed, and converge exponentially to the background state along distinct directions as curves in the 
distance- velocity plane, compare Figure fTTST c) and [looOOj . Near a turning point of the flux, solitons can 
travel with the sound speed evaluated at the background state. Such solitons converge algebraically 
to the background state as ^ ±oo and along the same line in the distance-velocity plane forming a 
cusp, see Figure HnT c). 

Concerning fronts it has been proven in [looOO] that fronts bifurcate from turning points r» of $' with 
$(4)(r,) ^ if and only if < q, i.e., convex-concave flux. These fronts travel faster than the 

sound speeds of left and right states, have monotone profiles and converge to the cndstatcs exponentially, 
compare also |HR09j . In our context fronts appear as conservative shocks. 



3 Numerical simulations of FPU-Riemann problems 

All simulations in this paper describe Riemann problems with cold initial data. That means for given 
left state wl = (fL, wl) and right state mr = (rR, dr) we initialise the atomic distances and velocities by 



(r„(0), «,(0)) = 



{r-L, vi) for ea < a^, 
(m, v-r) for sa > a* 



where e ~ 1/N is the scaling parameter and a* denotes the macroscopic position of the initial jump. We 
impose the boundary conditions 

VN+iit) = ujv(i), ro{t) = ri{t), 

so (|2.3p becomes a closed system for the 2N unknowns ri...rjv and vi-.-v^. These conditions are appro- 
priate since we start with piecewise constant initial data, and stop the simulation before any macroscopic 
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Figure 3.1: Numerical results for potential (|3.ip and Riemann initial data (|3.2p : Snapshots of atomic 
distances and velocity against the macroscopic particle index a for several macroscopic times and N = 
4000. On the macroscopic scale the solutions becomes self-similar: a left moving rarefaction wave and a 
right moving dispersive shock are separated by a cold intermediate state. 
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Figure 3.2: Atomic distances from Figure [XT] for different values of N . 



wave has reached the boundary. For the numerical integration of (|2.ip we use the Verlet-scheme, which 
is a symplectic and explicit integrator of second order [SYS97[ IHLW02j . The microscopic time step size 
At is independent of and small compared with the smallest inverse frequency of the linearised chain. 

3.1 Self-similar structure of solutions 

Typical examples for the numerical outcome of an atomistic Riemann problem are given in Figures 13.11 
and 13.21 where we plot snapshots of the atomic distances and velocities against the scaled particle index 
a = a/N. The used potential 

$(r) = exp (1 - r) ~ (1 - r) + ^{r - if (3.1) 

is a modified Toda-potential with strictly convex flux <i>' and the initial data are given by 

a, = 0.6, (rL, wl) = (0, 0), (rR, v^) = (0, 1). (3.2) 

In Figure [3TT] we fix N and plot the solution for i = 0, and t ~ 0.15, and t = 0.3, whereas Figure [321 
shows the numerical results for increasing N at the same macroscopic time t = 0.3. Recall that according 
to the hyperbolic scaling the microscopic time is always proportional to A^. 

The simulation indicate that the atomistic solutions indeed converge on the macroscopic scale to 
some limiting Young measure which is self-similar in t and a — a, . In the limit N —f oo we predict 
that the solution consist of a cold rarefaction wave and a dispersive shock which are separated by a cold 
intermediate state. In the cold regions the atomic data can be expected to converge to a macroscopic 
function, so in each point (Z, a) the Young measure is a delta distribution. In the dispersive shock, 
however, this measure is nontrivial but the envelopes (and likewise the local mean values) still converge 
to functions. 



Figures 13.11 and 13.21 provide of course a merely qualitative confirmation of our interpretation of the 
numerical data. A refined quantitative measurement with different A^ would be possible but requires 
much more numerical effort for the following two reasons. 
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(1) Since information propagates with infinite speed in the lattice system (|2.1[) cold states manifest 
only in the limit N —^ oo. For finite N we find small fluctuations everywhere due to the discreteness 
of a. Hcuristically, we expect the amplitude of the fluctuations to decay exponentially with outside 
the space-time cone spanned by the fastest macroscopic speeds, but only algebraically with TV inside 
this cone. This expectation is supported by rigorous results for the Toda chain and the harmonic chain, 
see |Kam93[ IMP09| , and implies that the cold intermediate state has superimposed fluctuations with 
amplitude 1/^/N. An accurate measurement of intermediate states and wave speeds inside the above 
cone therefore requires simulations with very large N. 

(2) Due to the oscillatory nature it is notoriously difficult to compare quantitatively the dispersive 
shocks for different values of or t: Accurate values for the evolution of the envelopes are very hard to 
measure numerically, and the precise values of averaged quantities such as local mean values or numerical 
distribution functions depend for finite N on the details of the implemented averaging algorithm. 

In this paper we focus on the qualitative properties of FPU Riemann solutions. We aim to understand 
the macroscopic selection and composition rules for elementary waves and how turning points of $ effect 
the qualitative structure of Riemann solutions. In particular, we do not intend to measure numerical 
convergence rates or to predict the wave parameters quantitatively. 

3.2 Elementary waves 

The following key observation about solutions to FPU-Riemann problems reflects the hyperbolic and 
modulation nature of the limit e ^ 0, but has not yet been proven rigorously. 

Observation 3.1. For strictly monotone and nonlinear flux $' where has at most one root in the 
range of the solution we observe the following. 

1. The macroscopic dynamics for cold Riemann data is self-similar and hence reducible to the macroscopic 
velocity variable c = (S — a*)/t. 

2. The arising measure at each [a, t) is either a point measure or supported on a closed curve that 
is generated by the distances and velocities of a wave train profile. Therefore, we can describe the 
macroscopic limit by a family of modulated wave trains parameterized by c. 

3. The macroscopic solution to each cold Riemann problem consist of a finite number of self-similar 
waves. These elementary waves are 

1. cold rarefaction waves, 

2. dispersive shock fans connecting two cold states, 

3. energy conserving jumps between two cold states (only for flux with turning point), 

4-. dispersive shock fans connecting a cold state and a constant wave train; these shocks always come 
as a counter-propagating pair. 

Structure of dispersive shocks Our basic numerical observations concerning dispersive shocks 
are as follows, see also Figure 11.31 As mentioned in |Jl] dispersive shocks have been studied for certain 
potentials and in other contexts, but we have not found the following explicitly mentioned for FPU. 

Observation 3.2. In the numerical simulation of cold Riemann problems dispersive shocks appear with 
oscillatory atomic data between two constant states, see Figure \1.3l Within a dispersive shock the self- 
similarity variable c = {a — oi^) /t ranges between the shock back velocity Cb and the shock front velocity 
Cf . The atomic oscillations have monotone envelopes with maximal amplitude at the front and vanishing 
at the back. There exist dispersive 1- shocks with Cf < and Cf < Cb as well as dispersive 2- shocks with 
Cf > and Cf > Cb. 

Moreover, our numerical results suggest the following fine-structure of the oscillations within a dis- 
persive shock. 

Observation 3.3. A dispersive shock consists of a one-parameter family of wave trains parameterised 
by c ~ {a — ol^)/t. Within the dispersive shock the parameter modulation is smooth (and hence follows a 
rarefaction wave of Whitham's modulation equations). Dispersive 2-shocks have the following properties 
(1-shocks accordingly due to symmetry). 
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Figure 3.3: Snapshots for the Toda chain with initial jump at a^, = 0.5 (see vertical lines), (a), (b) 
Distances and velocities for the classical 'supercritical shock problem' generating a permanent thermal- 
ization via binary oscillations; (c) Single dispersive shock where front and back counter-propagate; the 
asymptotic state at a* is a wave train with wave number ^ 0.47. 



1. The measure /i(cb) at the back of the shock reduces to the point measure generated by the constant left 
state Ml, and the local mean values of atomic distances and velocities smoothly connect to ml- 

2. The measure fJ.{c[) at the front of the shock converges to a soliton with background state mr, and the 
local mean values are continuous but not differentiable at Cf . 

3. The family of curves supp (/i(c)) with < c < Cf is nested. 

4- The dispersive shock is compressive in the sense that both positive characteristic speeds X± of the 
p-system point into the fan, i.e. A-|-(ml) > Cb and A+(ur,) < Cf . 

5. The Rankine-Hugeniot velocity of the jump lies strictly between Cb and Cf . 

For sufficiently small jump heights |rR — tl] + |wr, — ^lI the shock back and fr'ont move in the same 
direction, this means we have either Cf < Cb < or Cf > Cb > 0. The classical 'piston problem', however, 
shows that the situation is more subtle for large jump heights. The initial data in this problem describe 
an evenly spaced chain with positive left velocities and negative right velocities. It has been observed in 
[HS78j that sufRciently large ('supercritical') jumps in the velocity generate a transition from a pair of 
counter-propagating dispersive shocks with cold intermediate state to incomplete dispersive shocks whose 
'backs' are constant binary oscillations which replace the cold intermediate state. This phenomenon is 
related to dispersive shock fans with counter-propagating back and fr'ont. sec Figure [??^ c) . More precisely, 
passing the critical jump height from below the shock back velocities change their sign and the oscillatory 
intermediate state for supercritical data results from the interaction of two dispersive shocks. Analysing 
the back velocity of a dispersive shock might be a fruitful approach to determine the critical jump height 
in the piston problem. 



Conservative shocks in FPU chains As discussed in ^ numerical simulations indicate that 
conservative shocks appear in FPU Riemann problems only if they are supersonic. For illustration we 
consider the two quintic potentials 

3 4 5 

$(r -t- 2) = - — - — -f — (3.3) 
^ 6 24 120 ^ ' 

S 4 ^ 

lyt^ iv-i'-.' 



We plot the set of conservative shocks for these potentials in Figure [4?6l Note that due to Theorem l5.1l[ 5|) 
each of these sets consists of the diagonal and a closed curve crossing the diagonal at the turning points. 

The flux for (|3.3p has a turning point « 1.3 with $'^'(r») < (convex-concave), while the flux 
for (|3.4p has a turning point w 2.7 with $^^^(r*) > (concave-convex). Both potentials are convex in 
a neighborhood of r*, and the other turning points are outside the range of simulation. 

Potential (j3.3p allows for instance for the two supersonic conservative shocks 



rL = 2, rR«0.59, c^h = ±1.50, wr - wl « ±2.11, A±(rL) « ±1.41, A±(rR) « ±0.89. 
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+0.59 I 

0.85 1.0 0.8495 



Figure 3.4: Supersonic conservative 1- and 2-shock for potential 
indicated by the vertical lines. 



with N = 8000, t = 0.5 and a* 



In Figurc lB^ wc plot the solution to the corresponding FPU-Riemann problem (with ul = 0) and conclude 
that both the supersonic 1-shock and the supersonic 2-shock are captured by the atomic chain. 

The conservative shocks for potential (|3.4p that range over = 2.7, however, are subsonic. The FPU 
solution to 

rL = 4, rR, «1.24, c^h ~ ±1-46, i;r - « ±4.03, A± (rL) w ±1.83, A±(rR,) « ±1.73, (3.5) 

is plotted in Figure 11.41 and is far from a conservative shock. Recall that this is accordance with the 
non-bifurcation result for subsonic fronts in [looOO] . The solution in Figure 11.41 consists of a dispersive 
shock with attached rarefaction wave, that means both waves are not separated by a constant state. 
Consequently, the soliton at the front of the dispersive shock is no longer supersonic but sonic, i.e., 
its speed equals the sound speed of the background state. This is confirmed by the numerical data in 
Figure [LW c): the distribution function near the soliton has the predicted cusp shape. Compare with the 
non-degenerate exponentially decaying soliton in Figure fLST c) . 

4 Riemann solvers 

4.1 Towards an FPU-Riemann solver for the classical case 

In this section we describe an adaption of the classical p-system solver that accounts for dispersive shocks 
in the macroscopic solutions to FPU-Riemann problems. Based on Observations [3T2l and [3?3l we arrive at 
the following conjecture, which is illustrated Figure l4Tl 

Conjecture 4.1. From each state ul there emanate two dispersive shock curves I?_[itL] and 
with the following properties. 

1. Each state mr e 'D±[u]J\ can he connected with ul hy a single dispersive shock with sgn(cf) = ±1. 

2. The curves fit smoothly to the corresponding rarefaction curves TZ±[u]J\ (so that for small jump heights 
the dispersive and Lax-shock curves almost coincide). 

The macroscopic solution to FPU-Riemann problems with cold data and sufficiently small jump heights 
can he described hy the wave sets 

In particular, a FPU-Riemann solution consists of a unique 1-wave from W^^^[ul], an intermediate 
state um, and a unique 2-wave [um]- The intermediate state is cold if either a rarefaction wave 

occurs or if adjacent shock backs move away from each other; otherwise it is a wave train. 

To illustrate the first part of this conjecture we present simulations for the modified Toda potential 
(|3.ip . For given left state ul = (0, 0) and different values of tr < we choose wr such that ur = 
(''Ri ■f^R.) G '5-[itL], and study the macroscopic behaviour for the corresponding solutions to Newton's 
equations. The results for 5* ~ 0.5, t = 0.1, and = 4000 are depicted in Figure [4?2l For all values of 
tr we find a dispersive 1-shock whose front moves to the left, and an essentially cold intermediate state 
Um which gives a point in In all simulations there exists a right moving 2-wave, but this wave 

has much smaller amplitudes than the 1-wave. Hence, I?_ [ui] and iS_ [wl] are different but close to each 
other. Moreover, the simulations indicate the following behaviour for increasing jump height tl — fR. 
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Figure 4.1: Left: Sketch of the FPU Ricmann-solvcr for strictly convex (f>'. The wave sets yV^^^[uL] 
consist of rarefaction curves and dispersive shock curves, and decompose the plane into 4 regions DD. 
RD, RR, and RD. Right: The corresponding classical solver for the p-system with Lax shocks instead of 
dispersive shocks. 
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Figure 4.2: Three points from the dispersive shock curve I?_[(0, 0)] for potential (|3.ip . 
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Figure 4.3: Illustration of the FPU-Riemann solver for potential (jS.ip and initial data ()4.ip . The examples 
IVPl-3 correspond to the regions RR, DR, and DD, respectively, in Figure [4T] 



The shock front velocity ct decreases, whereas Cb increases and changes its sign at a critical value of tr. 
As mentioned, this critical value can be viewed as the analogue to the critical parameter in the piston 
problem, see Figure \3l3[ 

The resulting FPU-Riemann solver is illustrated in Figure 14.31 which shows the solutions to the 
following initial value problems with wl = (0, 0) and 

1: UR = (0, +1), 2: ur = (-1, 0), 3: ur = (0, -1). (4.1) 

It is important to note that the Lax shock wave sets S± [ul] generally differ from the dispersive shock 
wave sets Therefore, replacing by 2?±[ul] changes the Riemann solution, i.e., the precise 

values for the intermediate state and possibly the waves themselves. 

The difference between FPU chain and p-system can be quantified for the Toda potential (|2.2p . For 
the shock piston with ul = (1, +'2a), ur = (1, —2a), and < a < 1 (subcritical case) the classical solver 
for the p-system provides the intermediate state um = (tm, 0) with 

4a2 = {vM - l)(l - exp (1 - vm)) = {vm - if + 0{{rM - if) 
whereas the results in }Kam93j imply 

rM - 1 = 21n(a+ l)(a + 1)^ = 20 + 0(0^) 
for the FPU-Riemann solver. Both results are different but agree to leading order in a. 
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4.2 Riemann solvers for the p-system in the non-classical cases 



In this section we consider forces 3>' with a single turning point or where the solution ranges over at most 
one turning point. For the p-system, the classical Riemann solver cannot be used and also the above 
FPU-Riemann solver fails. To prepare the discussion of FPU-Riemann solvers in these cases, we first 
describe the relevant solvers of the p-system. 

The building blocks for each non-classical solver are modified wave sets which replace yy[uL] from the 
classical solver. Specifically, wave curves need to be adapted when intersecting the line r = r, , though 
these start out near wr = ul as classical wave sets (rarefaction waves and Lax shocks). The reason is the 
change in the sign of at r* , which implies that the Lax condition (jA.2|) or the compatibility condition 
A(rR) > A(rL) is violated. Note that the wave curves directed away from the turning point remain 
unchanged. For a single turning point of the flux the classical wave sets interact with the turning point 
as follows. 

Remark 4.2. Let ml with ri^ ^ be given and suppose convex-concave with $'^'(r<,) < 0. Then the 
curves 7?.+ [ul] and S-[ui] intersect the line r = r^ in the (r^v) -plane, whereas does not change its 
sign along 7?._[ul] <md S+[ui]. The same holds for concave-convex $' with $('^-'(r,) > if we replace + 
by -. 

Proof. Wc start with convex-concave $. For tl < we have $"'(rL) > and the formulas of Appendix 
\K\ implv that both 7?.+ [ul], iS_[ul] point into direction of increasing r, whereas r decreases along 7?._[wl] 
and iS+[wl]; compare Figure for an illustration of the wave sets of ul- The same is true for tl > r» 
as $"'(rL) < implies that now r decreases along 7?.+ [wl], 5- [ml]- Finally, the proof for concave-convex 
is analogous. ■ 




(a) (b) 



Figure 4.4: Modified wave sets of the conservative p-system solver: (a) modified shock curve, (b) modified 
rarefaction curve; C=compressive Lax shock, N=non-classical conservative shock, R=rarefaction wave; 
$'(r) as in the main text. In the insets we sketch an example for the different composite waves in each 
segment of the rR-axis. 

The conservative solver The modified shock curve 5f°"*'[wL] of the conservative p-system solver 
is illustrated in Figure I4.4r a) for convex-concave $' and tl > This is done by parameterizing 
^consjy^j by rpj and plotting $'(r) = <i>'(r) + p\r + p2 to give an illustrative graph. In accordance with 
Remark 14.21 the curve 5^°"'' [ml] starts out as for tr ^ tl- The wave set modification for large 

tl — r-R, goes via conservative shocks and works as follows. Along 5^°"''[ul] there exist a unique state 
^0 = (^Oj *^o) G 'H^[uy\ that can be reached from ml with a single conservative 1-shock. This state uq 
determines another state U2 = (^2, V2) corresponding to a Lax shock in 5_ [ul] such that at r2 the secant 
slope from tl to r2 coincides with the slope of the secant from tl to rp. This reads 

$'(rL) - $'(r2) $'(ro) - $'(r2) 

ro < r2 < TL and = , 

f L - f 2 r^- r2 

and implies that the Lax shock connecting ml to M2 and the conservative shock connecting ml to mq have 
the same Rankine-Hugeniot velocity. Note that this relation and in fact all conservative shock distance 
data are the same for and 
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The modified rarefaction curve 7?.™"''[rL] is illustrated in Figure H^ b) in the same way. The following 
lemma, which follows directly from results in [Smo941 ILeF02| , precisely describes the modified shock and 
rarefaction curves (for both convex-concave and concave-convex $'). 

Lemma 4.3. 

1. If^"'{r) 7^ for all r G Iq '■— (?'R,f'L),. then the solution coincides with the classical solution and 
consists of at most two uniquely chosen rarefaction fans or compressive shocks. 

2. Suppose $"'(r*) = for a unique and let wl = with tl > r^ be given. Then the following 
are uniquely defined for tr < tl- The solution ro of J{ri^,r^) ~ 0, the solution r,„ of J(r^,r„i) = 0, 
and the solutions ri of |crh(''L, ''i)! = |crh(?'L, ''m)| arid r2 of |crh(?'L, '"2)! = |crh(''L, ?'o)|- It holds that 
ri < To < r2 < < tl- Note that ro,r2 depend only on tl whereas r„i,ri are functions of r^ . 

(a) Let vq be such that uq := (?'o, Vo) G 'W±['"l]- ^ right state with tr < tl lies in the modified shock 
wave set 5™"^ [ml] i/wR G 5±[wl] for r2 < r^ , wr e 5±[iio] for ro < 7'r < 7-2, and ur e 7?.±[wo] for 
rn < rQ. By definition of ro , the shock from ul to uq is conservative, and the solution amplitude is 
discontinuous at r2. There always is an intermediate state between rarefaction fan or compressive 
shock and conservative shock. 

(b) Let Vm be such that u„i :~ (rrrnVm) G 'H±[u^]. A right state wr with tr < tl lies in the modified 
rarefaction wave set 7?.™"''[ul] «/ wr G 7?.±[ul] for tr > u,n € T^±[uh] for rg < tr < 

G S±[ul] for ri < tr < Tq, and mr £ H^[ui^] for tr < r^. 
By definition of r„i , for ri < tr < the shock from u„i to mr is conservative, and the solution am- 
plitude is discontinuous at ri . For tq < tr < r* the rarefaction fan is attached to the conservative 
shock, while for ri < tr < tq the compressive shock is not. 

3. If the flux $' has several turning points, then the modified wave sets are unchanged as long as uj^, ur 
are such that only one turning point lies in [ri,rm] and [rQ,r-if\. 

Proof. (1.) Lemma [5.3( 2.') shows that ^ in this case. Hence, the conservative Riemann solver 
coincides with the classical solver jLeF02| . For this solver, the p-system is solved uniquely in terms 
of at most two rarefaction or shock waves |Smo94j . (2.) The proof for the shock case is an immediate 
consequence of Theorem IV. 4. 3 (see also Theorem II. 4. 3) of |LeFd2, . The rarefaction case is not explicitly 
proven m |LeF02| . but follows from the exposition, cf. |LeF02| p. 163 (see also Theorem II.5.4). (3.) This 
follows from the independence of the solution on $ outside this range. ■ 

Due to this construction, the Riemann solution generically contains conservative shocks despite the fact 
that these are of higher codimension in the space of left and right states: A solution will consist of three 
elementary waves instead of two whenever a conservative shock is possible. Also note that the solution 
is non-monotone whenever a compressive and a conservative shock are connected in a solution. 

The dissipative solver We refer to the 'maximum entropy dissipation' solver in |LeF02| as the 
dissipative solver. This solver is much simpler than the conservative solver, and we do not explain it in as 
much detail. Wc plot the modified wave sets which determine the solver in Figure [4751 for concavo-convex 

and tl > r*. Compared with the conservative solver, the regions with conservative shocks have shrunk 
to points. The distance values where wave sets need to be modified are the turning point and the 
value Tg where the compressive shock has extremal velocity, i.e., where it coincides with a characteristic 
velocity. 

4.3 Towards an FPU-Riemann solver for the non-classical supersonic case 

We numerically tested the predictions of the conservative solver about the modified wave sets by simula- 
tions with initial data on the p-system shock and rarefaction curves. From the classical case we expect 
that compressive shocks are replaced by dispersive shocks in the FPU chain. Indeed, up to this modifi- 
cation, for convex-concave potentials the conservative solver qualitatively makes the correct predictions 
for FPU chains. 
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Figure 4.5: Modified wave sets of the dissipative p-system solver: (a) modified shock curve, (b) modified 
rarefaction curve; symbols as in Figure WM 
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Figure 4.6: Part of conservative shock data for the potentials (|3.3p in (a) and (|3.4|) in (b). The horizontal 
lines mark the turning point. Bullets mark the locations of (tx^^tyC) for the simulations the simulations 
plotted in Figures l4Jl and [48l for (a), as well as Figures [4?9l and [47TOl for (b). 



In order to illustrate the structure of the modified wave sets for (jS.Sp we proceed as follows. We fix 
the left state ul = (2, 0), and for the points Px marked in Figure B^ a) we solve two Riemann problems 
denoted by Sx and Rx. The value for is determined by Pec, whereas wr is chosen such that, for the 
p-system, 5*2; and Rx correspond to a single 1-shock and 2-rarefaction wave, respectively. The numerical 
results for the atomic chain are plotted in Figures 14.71 and 14.81 

Neglecting small waves and fluctuations (caused by the computational boundary and the positivity 
of e) the simulations provide numerical evidence that the macroscopic limit can indeed be described by 
a modified conservative solver. The solutions for 51 to 5*8 correspond to those in Figure I4.4f a) when 
replacing compressive by dispersive shocks, and we use the inset titles C, N+C, N+R in the following. 
For small jump, i.e. rR ^ tl, the chain produces dispersive 1-shocks with amplitudes proportional to 
the jump height (52, 53 = C), and there exists a critical point tr = f2 at which the unique supersonic 
conservative shock nucleates. For tr < f2 the conservative shock persists whereas the dispersive shock 
shrinks (54, 55, 56 = N+C) and transforms into a rarefaction wave (57 = N+R). Note that, in contrast 
to the conservative p-system solver, the nucleation of the conservative shock is a continuous transition in 
the envelope since the dispersive shock extends to the nucleating intermediate state. In the same way the 
solutions of R\ to R% transform to those of Figure H^ b). again using inset titles: i?l, i?3 = R, i?5, i?6 = 
R+N (with i?6 near N only), P7, i?8 ^ C+N. 

We summarize these modifications for the non-classical supersonic case in the following conjecture. 

Conjecture 4.4. For convex- concave flux Lemma \4.3\ holds under the following modifications and thereby 
defines the modified shock wave sets S^^[ui^] and the modified rarefaction wave sets TZjF^[uT_,] for the 
macroscopic FPU chain: (1) Replace Lax-shocks by dispersive shocks in the definition of S^^^^[u]J\, i.e., 
use 2?±[ul] from ' ^4-l\ (^) Replace r2 by the solution f2 of \c[{ri^,f2)\ = |crh('^Li '"o)!; md ri by the 
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Figure 4.7: Simulations for data on the supersonic 1-shock curve for potential p.3p and ml = (2,0), 
a* = 0.9, t — 0.5; N — 2000. The values for tr are those in Figure [TBT a). so 'IVP Sx' crosses the turning 
point for increasing x. 
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Figure 4.8: Simulations for data on the 2-rarefaction curve for potential (|3.3p analogous to Figure [47 



solution fi of |cf(r„,rL)| = |ci.h(ri, r™)|. 

4.4 Towards an FPU-Riemann solver for the non-classical subsonic case 

For subsonic conservative shock data, the situation is entirely different. Recall that in Figure [F4] subsonic 
conservative shock initial data did not produce a conservative shock. From ^ recall that fronts of 
the infinite chain do not bifurcate in the subsonic case. Indeed, simulations analogous to those in the 
supersonic case yield the drastically different results plotted in Figures 14.91 and 14.101 Since conservative 
shocks are absent and composite wave of rarefaction and dispersive shocks occur, we compare these with 
the solutions of the dissipative p-system solver. It is the only solver in the (natural) family of solvers 
studied in [LeF02| that does not use non-classical shocks, compare Figure H75l 

We investigate the modified wave sets for potential (|3.4p as before. For given left state ul = (4, 0) 
we choose several values of tr , see Figure and determine such that mr = (rR , wr ) belongs to 
7?._[ul] and The numerical results are plotted in Figures 14.91 and 14. 101 for comparison with the 

dissipative p-system solver we use the inset titles from Figure l4?5l For the shock initial data SI-S5 the 
chain produces single dispersive shocks {S1-S5 = C) with increasing amplitudes, decreasing back speeds 
Cb and increasing front speeds Cf. Here Cf is always larger than the corresponding Rankine-Hugeniot 
speed Crh and the speed of the conservative shock. For S6-S8 = C+K we find a qualitatively different 
solution with increasing rarefaction waves that are attached to the same dispersive shock. 

On the other hand, in the sequence R1-R8 the solutions are single rarefaction waves {R2 ^ R), from 
which a dispersive shock nucleates between i?3 and i?4 = R+N, when crossing the turning point r^, see 
Figure HTBl b). The rarefaction fan shrinks from R5 to R7, and eventually the solution consists of a single 
dispersive shock (not shown) corresponding to inset C in Figure B31 b). We conjecture that the solutions 
from both Figures and HTUI can be understood by modifying the dissipative solver as explained below. 
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Distances, IVP S1 Distances, iVP S3 Distances, IVP S5 




Figure 4.9: Simulations for data on the 2-shock curve for potential p.4p and ml = (4,0), and a* — 0.1, 
= 3000, t = 0.4. The values for tr, are those in Figure [4?6r b). so TVP 5a;' crosses the turning point 
for increasing x. The vertical lines indicate the locations of the 2-shocks of the p-system corresponding 
to (rL, tr), and the conservative shock corresponding to itL- 

Distances, IVP R2 Distances, IVP R3 Distances, IVP R4 




0. 0.5 0.9 1. 0. 0.5 0.9 1. 0. 0.5 0.9 1. 

Figure 4.10: Simulations for data on the 1-rarefaction curve for potential p.4p analogous to Figure SHI 



velocity 



nucleation 




To r2 ri^ ro r2 r* tl ^ 

(a) (b) 

Figure 4.11: (a) Sketch of the explanation for non- nucleation in the subsonic case. Here the front velocity 
Cf is larger than the nucleation velocity so that conservative shocks cannot be selected, (b) In the 
supersonic case, the nucleation cannot be missed in this way and occurs here at r2 ■ 



However, since there are no non-classical shocks to test against, and since we cannot compute Cf from 
the left and right states, our evidence is weaker than in the supersonic case. 

We numerically observed the absence of subsonic conservative shocks for various potentials. An 
explanation on the level of the Riemann solver is the following. Along the conservative shock-curve 
^^ons q£ ^-j^^ p-system, the nucleation of the conservative shock connecting to ul occurs when it is as 
slow as the compressive shock, i.e., Crh(''L,''o) = Crh(''L, ''r)- However, in the FPU case, the criterion is 
naturally modified to equality of conservative shock velocity and front velocity of the dispersive shock, 
i.e., Ci.h(rL,ro) = Cf(rL,rR). Recall that we have Cf(rL,rR) > Ci.h(rL,rR) for ?-r 7^ 7-l due to Observation 
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Figure 5.1: Sample sketches for J{ri^, tr) = 0. Between tl and the areas above and befow the secant 
hne up to the graph of <&' (shaded) are equal, (a) The secant transversely intersects the graph on left 
and right so that the set D is locally a monotone curve in the (rL, r-pj-plane. (b) Secant and graph are 
tangent at tl so that D has a local extremum in the (rL, 7'R)-plane. $' has (at least) two turning points 
in [rpijTL]. (c) The secant is tangent at both intersection points, and both tangencies point in the same 
direction, hence it is a local extremum of J and the point (rL,rR) is isolated in D. $' has (at least) four 
turning points in [rR,rL] 



I3.2i and note that both velocities converge to the characteristic velocity as — * tl- It is therefore 
plausible that ct(rL,7'p{) > Ci.h(?'L, tq) for all tr so that the nucleation criterion always fails. We sketch 
this situation in Figure [4. llT a). Note that the relative locations of the 2-shock and the dispersive shock 
front in Figure 14.91 support that this ordering indeed occurs for the potential (|3.4p . 

In contrast, in the supersonic case $('')(r,) < the velocity curve is unimodal with a maximum, 
so that nucleation cannot be missed in this way, see Figure I4.1ir b') . We thus arrive at the following 
conjecture. 

Conjecture 4.5. Let he the unique turning point of ^' . If ^^'^\r^) < 0, then non-classical shocks 
are absent and the solution is qualitatively according to the dissipative solver. More precisely, the wave 
sets depicted in Figure [775| need to be modified as follows: (1) Replace Lax-shocks by dispersive shocks in 
the definition o/5^™^[ul]; «-e., use 'D±\ui] from ^4-l\ f^j Replace r^ by the solution fg to Cf(rQ,rL)^ = 
$"(fo) and rl by the solution to Cf(fg,rL)^ = (f>"(rL). 

5 Properties of conservative shocks 

In this section we study conservative shocks, that means solutions to the three independent jump condi- 
tions (j2.7p . Eliminating the velocities ul, vr, an Cih one finds that each conservative shock is an element 
of 

D = {{ri^, tr) : Jin^, rR.) = 0}, (5.1) 

with J' as in (|2.8p . Conversely, each point in D defines both a conservative 1-shock and a 2-shock, which 
are unique up to Galilean transformation and differ only in sgn|t;]| = sgncrh- The geometric interpretation 
of J7 = is that the signed area between the graphs of and the secant line through $'(rL) and $'(rR) 
is zero over [fl , tr ] , compare Figure 15.11 

To characterise the structure of D in presence of several turning points of we use the notation 
CL |A±(rL)| and cr |A±(rR)|. 

Theorem 5.1. For $ G C''(M) the set D has the following properties. 

1. Off-diagonal data (tljTr) e D with rj^ ^ tr exists if and only if ^' has at least one turning point in 
the interval (7'l,?'r). 

2. Let I d R be any interval containing a single turning point of $' . Then D H I x I is the graph of a 
strictly decreasing function which crosses the diagonal {tl = tr} at the turning point. 

3. The conservative shocks corresponding to D C] I x I are undercompressive. If^^'^'^{r^) < 0, they are 
supersonic, and if <I>*^'*^ (r*) > subsonic. 
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Distances, IVP 1 Distances, IVP 3 




(a) (b) (c) 

Figure 5.2: (a) The set D for the potential (|5.2p . i.e., on solid curves holds J = 0: ofF-diagonal segments 
of D are labelled according to the ordering of c^j^ and \\. (b) Solution for a conservative shock with 
initial data at point (1) in (a), (c) Solution for the initial data from point (3) in (a); the solutions appears 
to contain a jump from wave train to wave train. In (b) & (c) the initial data ranges over two turning 
points of <!>'; N = 4000 initial jump = 0.5. 



4-. Compression changes precisely at local extrema in the coordinate directions of D in the (ri^^r^^) -plane. 
At extrema in the rj^-direction crosses Crh{rh,TYi)^ , and at extrema in the r^-direction c£ crosses 
Crh(rL, m)^- 

5. // $' has precisely two turning points, then D is the union of the diagonal {rL = tr} and a closed 
curve crossing the diagonal at the turning points. 

6. The set D does not have bounded connected components if $' has three or fewer turning points, and 
D \ {tl = rn} is bounded if the number of turning points is even. 

In Figure [STW a) part of the set D is plotted for the potential 

$(r + !) = ?' + -r^ + —r^ - - cos(2r) + — sin(3r). (5.2) 
2 20 4 10 

Changes in shock type occur for instance at point 2, which is an extremum in the r'R -direction so that 
becomes larger than cf.-^^ in the direction towards the nearest extremum in rL-direction at point 4. The 
conservative shocks on this curve start out supersonic, hence in the segment between points 2 and 4 the 
1-shocks are compressive and the 2-shocks are rarefaction shocks. At point 4 the term Cp becomes larger 
than c^jj when crossing it away from point 3, and so conservative shocks beyond point 4 are subsonic. 

Remark 5.2. 1. Conservative shocks do not have a preferred direction of propagation, and are isolated 
points on the shock curves S±[u-i\. In contrast, classical shocks have a selected direction of propagation 
in order to be compressive and generate continuous segments of S± [wl] • 

2. In lack of turning points, conservative shock data does not exist for the Toda potential and any cubic 
potential. For harmonic potentials, however, all shocks are conservative, and in fact contact discontinu- 
ities. 

3. For even potentials off-diagonal conservative jumps occur for tl = — ?'r,, because J{r, — r) = by 
symmetry. More generally, the symmetry $(r* + r) = <I>(r„ — r) with <I>"'(r*) = provides D <Z D with 
D = {(rL, ^r) = {r,: + r,r^ ^ r) : r G M} . Note that each conservative shock from D is degenerate as the 
characteristic velocities for left and right states equal. 

4. For quartic potentials (i.e., the classical FPU chains) all off-diagonal conservative data are given by D, 
so that non- degenerate conservative shocks occur only for potentials of polynomial degree five or higher. 



Proof of Theorem 15.11 
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1. This immediately follows from the geometric interpretation of J7 = 0, see Figure [01 

2. Again, the geometric interpretation shows that for fixed 7'l C / there is at most one solution r e / to 
i/('"L, f) = 0; similar for tr. Hence, the solution set is a monotone curve; that it decays follows similarly. 
Note that tangents of $ and the secant slope cannot coincide within /. 

3. & 4. These follow from Lemma [5.31 below. 

5. Consider (rL,r-R) as in Figure [STTJb) , where J{ri^,rii) = and the secant slope coincides with the 
tangent slope of at one end. Note that, since there are only two turning points, the graph of $' must 
lie on one side of the secant line near the tangency. Moving monotonically through the point where these 
slopes coincide, the area between the graphs can only vanish when the other point reverses its direction. 
Hence, the curve in the (rL, rR)-plane has an extremum, and this can only occur when both turning 
points lie in the interval (rR , tl). When continuing the curves from item 2, the tangency points of the 
graph intersections must be reached. Since there are no further changes in monotonicity, and curves are 
unique in suitable intervals /, the two curves emanating from the turning points must connect. 

It remains to show that there can be no bounded components of D that are isolated from the diagonal. 

Note that a stationary point (r£ , ) of requires tangency of the graph of $' with the secant segment 
at both r£ and rj^. Such a stationary point is a local extremum if the graph of $' is either above or 
below the common tangency at both r£ and r^, i.e., $"'(r£) and <i>"'(rp) have the same sign. For <I> with 
two turning points there exist a unique stationary point which is moreover a local extremum {r^,r^) 
with J^{r^,r^) 0, because the geometry implies that the enclosed area is only on one side of the 
secant. Hence, each zero of is a regular point and cannot be an isolated point of D. Now suppose for 
contradiction that a connected and bounded component of D existed. Then it must be a closed curve 
containing the local extremum (r£, rj^) in its interior. However, fixing r£ and moving from rfj towards 
r£ the secant segment stays above or below the graph of <J?' so that J ^ until tl = '"r on the diagonal, 
which is the contradiction. 

6. We continue the discussion of bounded isolated components from the previous item with a local 
extremum (r£, rj^) in the interior. Our arguments from above imply that the interval contains 
at least three turning points, and the tangency criterion for local cxtrema shows that the number of 
turning points in the interval must be even. 

Concerning boundedness of D\{rh ~ r^.}, note that for an even number of turning points, the convexity 
of outside a sufficiently large interval is the same. Hence, the secant line for sufficiently far distant 
tl and tr lies on one side of the graph of so that J ^ 0. 



Using the proof of the last item, it is not difficult to construct $ for which the set D consists of several 
bounded components that are disconnected from each other. The following lemma gives some more 
specific information. 

Lemma 5.3. 

1. Whenever ^'"{r^) = and $(^'(r*) 7^ for some 7% G M there exist a smooth locally unique curve 
ri^ !■ R(ri^) of solutions to ^ in {7'l 7^ tr,} U {(r*,r*)} and it has tangent (—1, 1) at (r,t,r*). 

2. Whenever c^[j(rL,7-R) 7^ Cr and J'(7'l,7'r) = then the set D from Theorem ] 5. 1\ is locally given by a 
function tr = R(ri^), which has a local extremum i/ c^h(7"Li ?'r) = c^- Similarly, whenever c^ij(7'Lj^r) 7^ 
Cl and Jiri^^r^) = then D is locally given by a function ri^ = R{r^), which has a local extremum if 
c?h('"L,?-R) = 4- 

3. On a curve (r, i?(r)) as in (1.) it holds that 

sgn(c^h - c±) = -sgn($''^'(r*)) for r « r*, sgn(cL - Cr) = +sgn(<I>'^^^ (r*)) for r « r*. 

Proof. 
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1. We readily compute that all first and second order partial derivatives of J vanish on the diagonal 
{tl = tr} and that all third order derivatives contain the factor ^'"(t'l/r,). Implicit differentiation of 
Jir^ R{r)) with respect to r then shows that bifurcations of solutions to = from the diagonal can 
only occur for $"'(r) = 0. The resulting bifurcation equation at r*. using x = i?'(r*), reads 

and has the solution x ~ I, corresponding to the trivial solution curve along the diagonal, and x ~ 
— 1, corresponding to the new bifurcating branch, as well as two complex conjugate roots that do not 
contribute to real solutions. 



2. Labelling variables as J^{ri,r2) = J^{ri^,r^.) we compute 



dr,J{n,r2) = —'^(rj)- 



Using this and the definitions of the velocities, implicit differentiation of J^{r, R{r)) = gives 

n/, N ^ dr,J{r,R{r)) ^ c,h(rL, tr)^ - cg^ 
^' dr,J{r,R{r)) c,h(rL,rR)2-4- 

The statement immediately follows from this formula. 

3. It follows from (1.) that i?(7'* + s) = — s to leading order, so that we can expand G±{s) 
c^jj(s* + s, — s) Cj_(r* =F s) and to leading order we obtain 

r* + s - (r, — s) 



2s ' ' ''22, 

Therefore sgn(c^jj — c\) = — sgn($(^' (r»)) for s w 0. The proof for — is completely analogous. 



The results of this section can be readily generalised to degenerate situations where ^'^"'^ = at the 
turning point: Theorem 15.11 is unchanged, but the bifurcation equations and local set of solutions in 
Lemma 15.31 does according to the degree of degeneracy. 

Finally, concerning FPU solutions for more than one turning point, consider the results for initial 
conservative 2-shock data that ranges over two turning points in Figure [??^ b) . f c) . The solution in (b) is 
the expected conservative shock, but the one in (c) is not. In fact, the solution in (c) appears to contain 
a (conservative) shock between wave trains which cannot be predicted from the p-system at all. 

Thus, the predictive power of the p-system investigated in this paper ends even qualitatively for 
more complicated fluxes. Extended systems such as the Whitham modulation equations arc required to 
understand the solution structures, but as even the hyperbolic nature of these is unknown, it is left for 
the future. 



A Classical Riemann solver for the p-system 

Formal substitution into (|2.4p of a self-similar ansatz in the variable c = (a — a*)/t gives 

— cf = i) , — ci) = $"(r)f, (A.l) 

where ' ~ d/dc, and this implies = $"(?■). The eigenvalues, i.e., characteristic velocities, and associated 
right eigenvectors of the p-system are given by A±(r) = ±-\/$"(r) and e±(r) = (1,A^), respectively. 
Strict convexity of (f> implies A_(r) < < \+{r) for all r, and the p-system is therefore (globally) strictly 
hyperbolic. Eigenvalues are genuinely nonlinear as long as docs not change sign. 

Concerning symmetries, the p-system respects Galilean transformations, and so the set of self-similar 
solutions is invariant under (r(c), v{c)) i— *■ (r(c), v{c) -\- vq) where vq is constant. Moreover, the p-system 
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exhibits the reflection symmetry that under c >—>■ —c each self-similar solution transforms according to 
(r(c), v{c)) 1-^ (r(-c), -w(-c)). 

The Lax theory for strictly hyperbolic systems with genuinely nonlinear eigenvalues is built up from 
the following two types of elementary waves. Rarefaction fans are smooth solutions, and (jA.ip implies the 
existence of two families, called 1-waves and 2-waves with c = A_(r(c)) and c = A+(r(c)), respectively. 
Shocks are discontinuous solutions, and satisfy (jA.ip only in the sense of distributions. This gives rise 
to the Rankine-Hugeniot conditions (|2.5p with shock speed Cih. The convexity of <f> and the identity 
'^rhlWI = ™ply the existence of 1-shocks with Crh < and 2-shocks with Crh > 0. Lax theory 

considers only compressive shocks that satisfy the Lax condition (with — and + for 1- and 2-shocks, 
respectively) 

A±(ml) > Crh > A±(ur,), (A.2) 

To describe the wave set yV[uL] for given left state wl = (fL, Vh), we define the integral curves ©-[wl] 
and C'_|_['Ul] by 



(r, v) e 0± [ml] ^ v^vlT V*"(s)ds, 

and the Hugeniot curves 7i_ [ul] and T-L+ [wl] given by 

(r, v) e 7i± [ul] ^ 1' = -fL T v/(*'(»') - ^'Ij^l))!?^ - tl). 

All these curves contain the point ml and can be globally parameterized by r. Due to (jA.ip we find 
-g^ = — c which implies 

(m, wr) e 7?-±[ul] <^ (r-R, wr) e Oii^L] and A±(rR) > A±(rL). 

The shock wave set = ^-[ul] U consists of all possible right states that can be 

connected with ul by a single Lax shock. By construction, any mr, must be an element of one of the 
Hugeniot curves of ml, but the Lax-condition (jA.2p selects one branch for each Hugeniot curve. 

Exchanging the role of left and right states provides sets 7?.[ur,] = 7?.-[mr] U 7^+[mr,] and 5[ur] = 
S- [mr] U S+ [mr ] , which contain all possible left states that can be connected to a prescribed right state 
MR by a single rarefaction wave or Lax shock, respectively. The standard Riemann solver is defined by 
these left and right wave sets as follows, see Figure HTT] For given (ul, ur) there exists a unique (see, 
e.g., |Smo94j ) intermediate state mm G W[ml] n W^[mr] such that mm S >V_[ul] and mr € >V+[mm]- The 
solution of the Riemann problem consists of the two elementary waves that connect ml to mm, and mm 
to Mr (one of these may be trivial). 

Shocks are classified as follows. A shock connecting ml to mr with speed Crh is called 

1. compressive, or Lax shock, if A(ml) > Cj-h > A(mr), 

2. rarefaction shock, if A(ml) < Cj-h < A(mr), 

3. supersonic, or fast undercompressive, if |ci-h| > |A(ml)| and jcrhl > |-^(mr)|, 

4. subsonic, or slow undercompressive, if |crh| < |-^(i*l)| and |crh| < |'^('i*r)|i 

5. sonic, if |crh| = |A(ml)| or Icrhl = |A(mr)|, 

with sgncih = sgnA < and sgncih = sgnA > for 1- and 2- shocks, respectively. All these definitions 
are invariant under reflections Crh <^ ^Cih, A(mr) <^ — A(ml). 
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